Atomic step motion during the dewetting of ultra-thin films 
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We report on three key processes involving atomic step motion during the dewetting of thin solid 
films: (i) the growth of an isolated island nucleated far from a hole, (ii) the spreading of a monolayer 
rim, and (iii) the zipping of a monolayer island along a straight dewetting front. Kinetic Monte Carlo 
results are in good agreement with simple analytical models assuming diffusion-limited dynamics. 



I. INTRODUCTION 

The shape changes of crystals are usually governed by 
the motion of atomic steps at the surface^* 2 -. The study 
of steps dynamics have permitted a better understand- 
ing of the morphological changes during various non- 
equilibrium processes, such as crystal growth 1 , surface 
electromigration 1 , or the decay of performed structures^. 
In the present paper, we focus on some key processes in- 
volving the dynamics of atomic steps during the dewet- 
ting of ultra-thin solid films. Dewetting is the process by 
which a continuous film breaks down into islands to lower 
the global energy. A large body of theoretical work has 
been devoted to the dynamics of dewetting with various 
continuum models^"—. Nevertheless, Kinetic Monte Carlo 
(KMC) simulation a 10 ! 11 have revealed that the dewet- 
ting process is controlled by the nucleation and dynam- 
ics of atomic steps on facets. Furthermore, recent studies 
have reported the observation of facets in dewetting ex- 
periments with Ag/Si(lll) 12 , SOI (Si/Si0 2 ) systemal&li, 
and YBaCuO films 15 , hence suggesting that atomic steps 
could play a crucial role in the dynamics. 

A full analysis of the morphological evolution dur- 
ing dewetting would require to investigate not only step 
motion, but also the formation of new steps via two- 
dimensional nucleation. Such an approach was proposed 
recently in Refs j 1Q i n . However, a quantitative descrip- 
tion of nucleation is in general delicate because it requires 
an accurate determination of the step free energy which 
enters in the nucleation barrier. We shall here discard 
the nucleation process, and focus on the motion of steps 
in various geometries. 

The paper is organized as follows. In section [III we pro- 
vide a description of the KMC model. Section Hill focuses 
on the layer- by- layer dewetting regime. In this regime, 
an island is first nucleated far from a pre-existing hole in 
the film. The island then grows by diffusion-limited mass 
transfer from the hole. We shall here present an analyti- 
cal solution of the island growth process, which provides 
a good agreement with the KMC simulations. The is- 
land finally grows so much that it surrounds the hole, 
and forms a monolayer 'rim' around it. In section HV1 we 
provide an expression for the dynamics of the monolayer 
rim once it has surrounded the hole. This expression is 
in quantitative agreement with the simulations, and con- 



tains some logarithmic corrections as compared to the 
simple linear behavior proposed in Ref 11 . In section W\ 
we consider a different regime where a thick facetted rim 
forms during the dewetting process. Islands nucleated 
on the rim facet then grow in size and spread along the 
dewetting front with a well defined velocity V Z i p . The 
numerical evaluation of this velocity is difficult because 
it would require delicate front-tracking procedures. We 
therefore present a method which allows one to deter- 
mine this velocity from a single snapshot of KMC simu- 
lations. This method, the results of which were presented 
in Ref. 11 , is reported here in details. 



II. KMC SIMULATIONS 

We model the dewetting of a crystalline film using a 
solid-on-solid model on a 2D square lattice. The lat- 
tice unit is a. The substrate surface is perfectly flat and 
frozen. The local height is z. On the epilayer atoms 
z > 0, and on the substrate z = 0. 

We employ KMC simulations to implement the dy- 
namics. Epilayer atoms hop to nearest neighbor sites 
with rates r n when they are in contact with the substrate 
(z = 1), and v n when they are not in direct contact with 
the substrate {z > 1). In our model, an atom needs to 
break all its bonds to hop. The hopping barrier is there- 
fore given by the binding energy of the atom. Hence: 
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where v$ is an attempt frequency, T is the temperature 
(in units with fog = 1), n is the number of in-plane near- 
est neighbors of the atom before the hop, J is the bond 
energy, and Es is the adsorbate-substrate excess energy. 
The model is presented in Fi^JTJ Note that we do not 
write explicitly the bond J between the moving atom 
and the atom directly underneath it in the energy bar- 
rier for atom motion in Eqs. (jll2p . Indeed, since this bond 
is present below any atom before the hopping event, we 
rather include its contribution exp[— J/T] in the prefac- 
tor vq. We choose J as the energy unit, so that J = 1 in 
the following. 

The algorithm used in the simulations is the following. 
We list all atoms into classes. Each class is characterized 
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by the number of in-plane neighbors n of the atom and 
by the existence or not of a nearest neighbor belonging 
to the substrate. At a given time t, we calculate the 
probabilities per unit time W{ of all possible events ( an 
event is the motion of an atom originally at position i 
) given either by Eq.(pQ) or Eq.(j2j), and the sum W of 
all those rates (for all mobile atoms). We increment the 
time by a 5t, which is equal to the inverse of the sum 
of the rates of all possible events, 1/W. This choice for 
St corresponds to the average value of the waiting time 
between two successive events (Poisson processes)^. We 
choose the event with probability Wi/W. The chosen 
atom moves with equal probability in any of the four 
possible directions. 

The parameter Es controls the wetting properties of 
the film on the substrate. In order to relate more 
precisely Es to surface energetics, let us perform two 
Gedankenexperiments. The first one consist in splitting 
an epilayer crystal into two parts perpendicular to the 
z axis. We obtain two surfaces, and the energy balance 
is: J = 2Eavi where Eav is the energy per site of the 
created facet. The second Gedankenexperiment is to sep- 
arate the substrate and the epilayer. The energy balance 
now reads: Eas + J — Es = Eav + Esv, where Eas 
and Esv is the epilayer-substrate and substrate- vacuum 
energies per lattice site. Combining these two relations, 
we find 



Es = Eav + Eas — Esv- 



(3) 



Since the simulations were performed at a fixed tempera- 
ture T, we should in principle consider a balance between 
free energies, rather than energies. However, since the 
substrate is frozen, the epilayer-substrate, and substrate- 
vacuum interfaces do not exhibit any configurational en- 
tropy. Their energy is therefore equal to their free energy. 
In addition we will work below the roughening tempera- 
ture, so that the epilayer- vacuum free energy can also be 
approximated by an energy. As a consequence, the pa- 
rameter Es is, to a good approximation, a balance of free 
energies. This parameter is then identical (with opposite 
sign) to the thermodynamic work of adhesion S defined 
in Ref JLia. 

When Es < 0, we have complete wetting. The regime 
of partial wetting is obtained for Es > 0. At low tem- 
peratures, the equilibrium island shape has a square base 
of lateral size L eq , and a height h eq , with h eq /L eq — Es. 
When Es ^> 1, the energy minimization will favor high 
islands. As Es is decreased to lower values, the equilib- 
rium aspect ratio of the island decreases, and the islands 
become flat for Es <C 1. 




FIG. 1: Kinetic Monte Carlo (KMC) model, with hopping 
rates r n and v n . 



mimics dewetting initiated at heterogeneous nucleation 
sites. When Es is small or h is large, dewetting occurs 
in a layer-by-layer fashion. The precise criterion for this 
regime to occur was given in Ref. 1 — . The process starts 
with the nucleation of a monolayer island on the film. 
Since this nucleation event is randomly located 11 , the is- 
land hole distance is usually large, as shown in Figj2jb). 
The island and the hole then grow simultaneously due 
to mass transfer from the hole to the island. We shall 
here focus on this island and hole growth process, and 
propose an analytical expression for the dynamics. We 
assume that the hole-island pair is isolated on the sur- 
face: no other hole or island exists in its neighborhood. 
In the initial stages, the distance between the hole and 
the island is larger than the typical sizes of the hole and 
the island. In this limit, we may assume that the hole 
and the island are circular. We also assume that the 
dynamics is diffusion-limited, as in RefiA!. From the as- 
sumption of fast adatom attachment-detachment kinetics 
at the island and hole edges, the concentrations is fixed 
to its equilibrium value. Hence, the concentrations is c^ q 
at the hole edge, and cf q at the island edge. In such a 
geometry, the diffusion problem can be solved in bi-polar 
coordinates 19 . The adatom mass flux from the hole to 
the island is then found to be: 



J 



2ttDAc 



\n[ird 2 /(AB) 



1/21 



(4) 



where Ac = c{ 



c^, and A and B are the areas of the 



hole and of the island, and d is the distance between the 
center of the hole and the center of the island. From mass 
conservation, the flux J is then related to the evolution 
of the hole area: 



_ ft 27rftDAc 
* ~h ~ hln[dH/(ABy/2y 



(5) 



where ft = a 2 is the atomic area. The hole and island 
areas are then related via global mass conservation 



III. ISLAND AND HOLE GEOMETRY 



d t B = hd f A. 



(6) 



A. Model 

Let us consider the dewetting process starting from a 
pre-existing hole in a film of thickness h. Such a process 



Since the island forms after the hole, we may assume that 
at t = 0, A = Aq and B = 0, so that 



B 

A= -r + A . 
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(7) 
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(b) 
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FIG. 2: Typical dynamics in the layer- by- layer dewetting 
regime. Parameters: T = 0.4, Es — 0.25, system size 
400 x 400, and h — 3. (a) A hole is performed artificially 
in the film, (b) An island nucleates far from the hole, (c) Dif- 
fusion limited mass transfer from the hole to the island leads 
to the simultaneous growth of the island and the hole, (d) 
The monolayer island surrounds the hole, (e) The monolayer 
island closes around the hole and forms a monolayer rim. (f) 
The monolayer rim grows and becomes circular. 



Using this relation into Eq.(j5j), we find: 

2ttQDAc 



d t A 



(8) 



h\n[d 2 7r/(hA(A- Ao)) 1 / 2 } 
the solution of which can be written in the implicit form: 
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(9) 



where A = A/ird 2 , and Ao = Ao/nd 2 . At long times, but 
when the hole is still small as compared to the hole-island 
distance, i.e. when 1 ^> A ^> Ao, one has: 
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B. Comparison to KMC simulations 

In order to compare Eq.(j9]) with KMC simulations, we 
need to evaluate the parameters which enter in the model. 
We choose the lattice spacing as our unit length, and 
the time unit is the inverse of the attempt frequency. 
Hence ft = 1 and vo = 1. Since an atom can move 
randomly to the 4 nearest neighbor sites, the diffusion 
constant on the film is D = 1/4. Moreover, we obtain 
the equilibrium concentration on the film c eq = e -2J / T 
from the detailed balance of attachment-detachment at 
kink sites. Since the curvature of the island is larger than 
the curvature of the critical island size, we may assume 
that the Gibbs-Thomson contributions proportional to 
the curvature, are negligible, and cf q « c eq . 

The chemical potential in the vicinity of the film edge 
is Es/Th, as shown in Appendix [A] Our reference sate 
is the state of equilibrium on top of the film, with equi- 
librium concentration c eq (this is the equilibrium concen- 
tration in the vicinity of an atomic step on the top of the 
film). Using the standard formula: 



Tln^, 



we find that 



^eq ' 



^Es/Th 



We therefore have 



Ac = 



-2J/T n E s /(hT) 



■)■ 



(ii) 



(12) 



(13) 



Using these parameters, the solution of Eq.(|9]) is in good 
agreement with KMC simulations, as shown in Fig|3l 
Note that the Eq.(j9]) cannot lead to a perfect agreement 
with KMC simulations due to inherent assumptions of 
the model. Indeed, we have assumed that the island and 
hole are isolated, neglecting interaction with the periodic 
images of the island and the hole. In addition, the con- 
dition that the island-hole distance d is much larger that 
the radii of the hole and the island is verified at the be- 
ginning, but not at end of the KMC simulation. Finally, 
we observe in Fig J3] that the condition A ^> Ao is not ver- 
ified in the KMC simulations. Therefore, the asymptotic 
linear solution ([TO]) cannot provide a good approximation 
for the KMC simulations result. 



IV. GROWTH OF A MONOLAYER RIM 
AROUND A HOLE 

A. Model 



Due to the smaller thickness of the monolayer as com- 
(10) pared to that of the hole, the island area increases faster 
than that of the hole, as shown by Eq.((6]). Once the is- 
land becomes large enough, it collides with the hole and 
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FIG. 3: Dynamics of one hole exchanging mass with a distant 
monolayer island. The symbols are the results of KMC simu- 
lations with: T,E s ,h = 0.4, 0.2, 2 (squares), and 0.35, 0.25, 2 
(circles), in a 400 x 400 periodic lattice. (In order to fit both 
curves in the same plot, the timescale of the first data set 
(squares) was multiplied by 6). The solid lines represents the 
solution of Eq.j9]) where the initial conditions and the value 
of d were measured from the simulations. 



surrounds it. Then, a different regime can be observed, 
where the hole is surrounded by a monolayer rim. This 
monolayer rim is initially very far from being circular. 
However, the diffusion limited transport provides more 
mass to the parts of the monolayer rim which are closer to 
the hole, resulting in a stabilization of the rim, which be- 
comes approximately circular. In this section, we present 
the exact solution of the circular rim evolution dynam- 
ics, and compare this solution to the asymptotic solution 
given in Refii. 

Mass exchange occurs on top of the monolayer rim via 
adatom diffusion. The adatom concentration field on the 
rim obeys: 



-d r \rd r c] — 0, 
r 

which has a solution of the form 
c = c In [r/r ]. 



(14) 



(15) 



We once again assume instantaneous attachment- 
detachment kinetics, so that c = c\ at r = R\ at the 
edge of the hole and c = c 2 at r = R 2 at the other edge. 
Moreover, we must impose mass conservation 



d t R 2 = ~^Dd r c\ 2 , 

d t Ri = -—Ddrcl^ 
rii 



(16) 



where hi = h + 1 is the rim height, and h 2 = 1 is the 
thickness of the monolayer rim. Introducing the areas 
A\ = irRl and A 2 = ttR 2 , we find 



d t Ai 



AttDQAc 
hi\n[A 2 /Ai 



(17) 



where Ac = c\ — c 2 is given by Eq. (fl3|) with the substi- 
tution h —J- hi. The evolution of A 2 is fixed by the mass 
conservation relation 



which implies 



h 2 d t A 2 = h 1 d t A 1 , 
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where A® are the areas at t = 0. Substituting this ex- 
pression for A 2 into Eq. ([T7]h one finds the solution in an 
implicit form: 

h 2 (A 2 \n[A 2 ]-A° 2 \n[A° 2 ]) 
-hi(Ai ln[Ai]-A?ln[A?]) 
= ATrDQAct. (20) 
At large times when Ai ^> A\ and A 2 ^> A 2l we find 



Aihi In 



w AttDQAc t. 

(21) 

Keeping the dominant terms only, we obtain a linear be- 
havior at large times 

AttDQAc 

Al w hiHhi/h^ (22) 

which is the asymptotic form presented in Ref. 11 . 



B. Comparison to KMC simulations 

In Fi^JU we have compared the KMC simulations re- 
sults with Eq. (|2Q|) and Eq. (|22|) . The agreement between 
the linear asymptotic formula ([22]) and KMC simulations 
results requires a fitting parameter (e.g. the reference 
time for which A\ extrapolates to zero). However, due 
to the divergence of the subdominant terms in Eq. ([2T]h 
the meaning of this parameter is unclear. 

Therefore, it is more satisfactory to use the full so- 
lution ([20]) . Indeed, the full solution requires no fitting 
parameter: we simply need to measure A\ and A 2 at an 
initial time once the monolayer rim is formed. The full 
solution is shown to be in good agreement with the KMC 
simulations in FigHJ 



V. ZIPPING OF A MONOLAYER ON A 
FACETTED RIM 

We now turn to another regime, which appears for 
large Es and small h. In this regime, a thick multilayer 
rim forms. The top of this rims is facetted. The increase 
of the rim height proceeds via the nucleation of mono- 
layer islands on the rim facet. Once they are nucleated, 
these islands grow in size so as to invade the whole rim 
facet, as shown on Figj5] Here, we want to focus on the 
growth process of the monolayer islands. 
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FIG. 4: Diffusion limited regime. Dynamics of A\ for a hole 
with a monolayer rim. The dashed lines correspond to the 
asymptotic linear formula Eq. (|22p . The solid lines correspond 
to the full solution Ea. (|20)) without fitting parameter. 



A. Monolayer zipping velocity 

Let us consider a straight dewetting front of height 
hi, with a rim facet of width t = X2 — xi. We aim to 
analyze the velocity of an island zipping along the dewet- 
ting front. We assume that the dewetting front does not 
deform and remains straight during this process. The 
driving force for the atoms to detach from the front, and 
to attach to the monolayer step is the difference between 
the chemical potential of the straight film edge E$/hi 
(see Appendix A), and the chemical potential of the step 
Qjk, leading to: 



A// : 



hi 



(23) 



where k is the step curvature, and 7 is the step stiffness. 
Here it is assumed that the temperature is high enough 
so that steps properties are isotropic. As a consequence, 
we have 7^7. 

The mobility of the atoms for diffusion between the 
film edge and the atomic step is M = Dc eq / (AT), where A 
is the typical distance between the fronts. The only avail- 
able geometric scale for a monolayer rim zipping along a 
dewetting front is the step curvature, and we therefore 
expect k ~ 1/A. 

As a summary, the zipping velocity should take the 
form: 



MA /1 



nDc e 



hi 



(24) 



We assume that the selected curvature is the one which 
maximizes the velocity: 



E s 
2hiQj' 



(25) 




x l x 2 

FIG. 5: Zipping of a monolayer island on the rim facet. Detail 
of a KMC simulation on a 1000 x 1000 lattice with T = 0.4, 
Es — 0.5, h = 3. On this image, hi — 8. 



Thus, the zipping velocity takes the form 

Dc ea E% 



V Z7 



c z 



Thh ' 



(26) 



where C Z i P is an unknown number. 



B. Time for the top step to cross the rim facet 

Two-dimensional islands are usually nucleated close to 
the front h\. The islands grow in the y direction with ve- 
locity V Z ip, as discussed in the previous subsection. But 
islands also grow in the x direction. In order to model the 
growth in the x direction we use a ID model, where a step 
is parallel to the front. This ID model only makes sense 
in the limit of small slopes d y x <C 1. We assume that 
dynamics is limited by diffusion on the top of the mono- 
layer island. The chemical potential difference between 
the dewetting front and the island edge is Es/(hi + 1). 
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Assuming once again fast attachment- detachment kinet- 
ics, the concentration at the step edge and at the front 
edge are respectively c eq , and c eq exp[Es/(hi + l)T\. The 
concentration profile is therefore linear, and the flux per 
unit length from the front to the island edge is: 



J 



Da 



11 ( e £?s/(/n + l)T _ 1 



)■ 



(27) 



where x s is the distance between the dewetting front and 
the island edge. The velocity of the island edge in the x 
direction is therefore: dx s /dt = QJ. From this equation, 
we find 



2(lDc eq (eEs/(h 1+ i)T _ ^ t 



1/2 



(28) 



The step will reach the step bunch at x s = £ in a time 
tg. For weak driving forces, E$ <C Thi, one has 



(29) 



C. Zipping length 

The direct observation of V z i p or tg in the KMC sim- 
ulations would require delicate front tracking procedures 
and we shall look for a quantity which can be observed 
on a single snapshot. From relation ([29]) . the tip length 
of the zipping monolayer reads: 



L zip 



E S 2 hi + 1 
zip 2^ hi ■ 



(30) 



The length L Z i P is the length between the tip of the zip- 
ping monolayer (in the vicinity of the front of height 
hi), and the point where the step edge of this monolayer 
reaches the bunch of height /12 (the distance is measured 
along y). 

Since the nucleation barrier for the formation of new 
monolayer islands grows with hi, the typical distance be- 
tween island nucleation sites increases exponentially with 
hi, as shown in Ref*ii. Thus, the length L Z i P is always 
smaller than the typical distance between island nucle- 
ation sites. Therefore, for large enough hi, the islands 
first reach the total width £ of the facet, and then zip 
along the dewetting front with a tip length L Z i P . This is 
indeed the scenario observed in KMC simulations. 



D. Comparison to KMC 

We have performed KMC simulations, with Es = 0.5, 
T = 0.4, h = 3, in a 1000 x 1000 system with total 
simulation time t = 3.5 x 10 8 . The simulation time is 
chosen to be large enough for the system to reach the 
late stage regime where L Z i P ^> £. Measurements were 
performed on five KMC simulations snapshots, each time 
measuring L Z i P from an average of its value at the two 



hi 


£ 


L z ip 


C zip 


6 


56 


86 


0.24 


6 


56 


54 


0.15 


7 


69 


150 


0.31 


7 


85 


180 


0.24 


8 


83 


124 


0.20 



TABLE I: Observed values of £ and L Z i P from KMC simula- 
tions with T = 0.4, E s = 0.5, and h = 3 in a 1000 x 1000 
lattice. 



ends of the expanding island on the top facet. From each 
measurement, we have extracted the value of C Z i P from 
Eq.([30j) with 7 = 0.43lfiJJ-2tt. The results are reported 
in Table I. We obtain the average value 



C zip « 0.25 ± 0.05. 



(31) 



The error here provides an indication of the order of mag- 
nitude for the observed variations of C Z i P extracted from 



different measurements. The result (|3T|) does not include 
the case where hi = 6, for which the condition £ <C L Z i P 
is not verified. Nevertheless, we see in Table U that the 
results with hi = 6 would provide similar values for C Z i p . 

Finally, we would like to point out that we have not 
yet studied systematically all possible regimes of zipping, 
due to limitations in the our KMC simulations. Indeed, 
for smaller values of Es, the simulation time is too large, 
and for larger values of Es, the solid on solid restriction 
(which forbids overhangs) should be released. Important 
improvements in our numerical approach are therefore 
needed in order to address these questions. 



VI. CONCLUSION 

In conclusion, we have analyzed the motion of atomic 
steps during the dewetting of ultra-thin solid films with 
various geometries: (i) simultaneous growth of a mono- 
layer island with a hole; (ii) expansion of a hole with a 
monolayer rim; (iii) zipping of a monolayer island on the 
top facet of a thick rim. 

Our KMC results are in good agreement with step 
models where mass transport is limited by diffusion on 
terraces between steps. These results demonstrate that 
the step model approach provides a systematic way to 
investigate the evolution dynamics of thin films during 
the dewetting process. 

Several additional ingredients could be added to the 
step models, such as anisotropy, non-trivial attachment- 
detachment kinetics (e.g. an Ehrlich-Schwoebel barrier, 
or step transparency), or elastic interactions between 
steps. We hope to report along these lines in the future. 

We acknowledge support from " nanomorphogenese" 
and "DeFiS" ANR-PNANO grants. 
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Appendix A: Local chemical potential at the edge of 
a facetted film 



Our variable here is xo(y), the position of the edge. We 
therefore have: 



The total free energy of a film of height ho at x > xo (y) 
(the height is zero at x < xo(y)) is: 



ST 



5x (y) \Sx (y) 



5AT 



T = J dso j e dge - Es& 1 J dy x (y), (Al) Using the above expressions: 



where so is the arclength along the edge of the film, and 
ledge is the free energy per unit length of the edge of the 
film. The total number of atoms in the film is: 



(A2) 



N = M -n 1 h j dyx (y), 

where A/o is a constant. Equilibrium reads: 
5 (J - - (jN) = 0. 



(A3) 



Eg 

li ^70^0 + -7—, 



(A4) 



(A5) 



where 70 = (ledge + ledge) / h o is tne stiffness per layer 
of the film edge (If the edge is a bunch of steps, and if 
the step-step interaction is negligible, then 70 is the step 
stiffness). Note that this derivation of the chemical po- 
tential relies on the postulate that the film profile across 
the film edge does not vary along along the film edge. 
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